#### This script will overlay response curves of species to microclim and bioclim
#### Need to include variable importance on each plot


### Establish directories

base.dir = '/home1/99/jc152199/MAXENT/output/'
out.dir = '//home1/99/jc152199/ChapterOne/Figures/ResponseCurves/'

### Read in raw physiology data from base.dir

minmax = read.csv('/home1/99/jc152199/MAXENT/output/PhysData/minmax.csv', header=T)

### List of species

species = list.files(base.dir)[c(1:length(list.files(base.dir))-1)]
species = species[-7]

### Loop through plotting for each species

for (spp in species) # Loop through species, leaving out COPORNA for now because of split between AU/SU

	{
	
	### Subset minmax data for one species
	
	spp.mm = minmax[which(minmax$Species==spp),]
	
	### List of all .dat files 
	
	mc.files = list.files(paste(base.dir,spp,'/microCLIM/output/plots/',sep=''),pattern='.dat',full.names=T) # List of all files in microclimate plot directory NB need to remove even file numbers
	bc.files = list.files(paste(base.dir,spp,'/BIOCLIM/output/plots/',sep=''),pattern='.dat',full.names=T)
	
	### Subset to just the 'only' .dat files
	
	mc.files = mc.files[grep('_only',mc.files)] # Data for plotting response curves
	bc.files = bc.files[grep('_only',bc.files)]
	
	### Read in summary data with variable importance
	
	mc.sum = read.csv(paste(base.dir,spp,'/microCLIM/output/maxentResults.csv',sep=''), header=T)
	bc.sum = read.csv(paste(base.dir,spp,'/BIOCLIM/output/maxentResults.csv',sep=''), header=T)
	
	### Open the png device driver
	
	png(paste(out.dir,spp,'_ResponseCurves.png',sep=''),units='cm',width=18,height=18,res=500)
	
	### Configure plot space
	
	par(mfrow=c(1,1))	

	### Position tracker
	
	h=1
	
		for (varx in c('05_')) # Loop through response for each surface, include underscore so that search for '15' doesn't find every file
		
			{
			
			### Remove underscore for naming
			
			var.names = gsub('_','',c('05_'))
			
			### Select the appropriate .dat files for the variable of interest
			
			mc.file = mc.files[grep(paste(varx),mc.files)]
			bc.file = bc.files[grep(paste(varx),bc.files)]
			
			### Get out contributions for the variable of interest
			
			mc.contrib = round(mc.sum[grep(paste(gsub('_','',varx)),names(mc.sum))[1]],1)
			bc.contrib = round(bc.sum[grep(paste(gsub('_','',varx)),names(bc.sum))[1]],1)
			
			### Read in the .dat files
			
			mc.data = read.csv(paste(mc.file,sep=''))
			bc.data = read.csv(paste(bc.file,sep=''))
			
			#### Establish limits based on bc data and mc data
	
			x.lims = range(range(mc.data$x),range(bc.data$x))
			x.lims = c(x.lims[1]-1,x.lims[2]+1)
			x.lims = round(x.lims,0)
			
			### Calculate some polygons that represent CTMin and CTMax 95% CI's (based on normally distributed data)
			
			if(nrow(na.omit(spp.mm))!=0)
			
			{
			
			max.poly = data.frame(x=c(spp.mm[1,5]-spp.mm[1,6]*1.96,spp.mm[1,5]-spp.mm[1,6]*1.96,spp.mm[1,5]+spp.mm[1,6]*1.96,spp.mm[1,5]+spp.mm[1,6]*1.96),y=c(1,0,0,1))
			
			### Establish limits based on bc data and mc data
	
			x.lims = range(range(mc.data$x),range(bc.data$x), range(max.poly[,1]))
			x.lims = c(x.lims[1]-1,x.lims[2]+1)
			x.lims = round(x.lims,0)
			
			}
			
			### Establish plot space and plot mc data as lines
			
			plot(mc.data$x,mc.data$y,xlim=c(x.lims[1],x.lims[2]),ylim = c(0,1) ,ylab = 'Probability of Presence',xlab = 'Mean Max Temp of Warmest Period',main=paste(spp,sep=''),type='l',col='blue') # Set the plot space and plot the first density object
			
			### Plot bc data as lines
			
			points(bc.data$x,bc.data$y,col='red',type='l')
			
			### Add a legend with variable contribution included
			
			legend('topleft',legend=c(paste('Summarized BRT Contribution ',mc.contrib,sep=''),paste('BIOCLIM Contribution ',bc.contrib,sep=''),'CTMax 95% CI'), text.col = c('blue','red','black'), lty = c(1,1,1), lwd=c(2,2,8), col=c('blue','red','#FF000025'), cex = .8, bty='n')
			
			
			### Extract n-value from phys data
			
			n.min = spp.mm[1,4]
			n.max = spp.mm[1,7]
			
			if(nrow(na.omit(spp.mm))!=0) # Conditional plotting
			
			{
			
			### Plot the polygons
			
			#polygon(x=min.poly[,1],y=min.poly[,2], col='#0000FF25', border=NA)
			polygon(x=max.poly[,1],y=max.poly[,2], col='#FF000025', border=NA)
			
			### Plot some text labels
			
			#text(x=min.poly[1,1],y=min.poly[1,2],label=paste(n.min))
			#text(x=max.poly[1,1],y=max.poly[1,2],label=paste(n.max))
			
			}

			### Plot lines to show the mean CTmin and CTmax
			
			#points(c(spp.mm[1,2],spp.mm[1,2]),c(0,1),type='l',lty=2, col='#0000FF50')
			points(c(spp.mm[1,5],spp.mm[1,5]),c(0,1),type='l',lty=2, col='#FF000050')
			
			### Adjust the position tracker
			
			h=h+1
			
			}
			
			dev.off()
			
			cat('\n',spp,' Image Completed\n')
	
	}
	
# Close loop
